Modelling power-law distributed interevent times 
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Many human-related activities show power-law decaying interevent time distribution with exponents usually 
varying between 1 and 2. We study a simple task-queuing model which can reproduce this property and we 
show exact results for the asymptotic behaviour of the model. The model satisfies a scaling law between the 
exponents of interevent time distribution (j6) and autocorrelation function (a): a + p = 2. This law is general 
for renewal processes with power-law decaying interevent time distribution. We conclude that slowly decaying 
autocorrelation function indicates long-range dependency only if the scaling law is violated. 

PACS numbers: 



Studying human activity patterns is of central interest due to 
the wide practical usage. Understanding the dynamics under- 
lying the timing of various human activities - such as starting 
a phone call, sending an e-mail or visiting a web-site - are 
crucial to modelling the spreading of information or viruses 
|Q]]. Modelling human dynamics is also important in resource 
allocation. It has been shown that for many human activi- 
ties the interevent time distribution follows a power law with 
exponents usually varying between 1 and 2. Processes with 
power-law decaying interevent time distribution look very dif- 
ferent from the Poisson process which has been used to model 
the timing of human activities before |Hl[3|. While time series 
from the latter look rather homogeneous, the former processes 
produce bursts of rapidly occurring events which are separated 
by long inactive periods. 

Some examples where power-law decaying interevent time 
distribution has been observed are email-communication 
(with exponent /3 ss 1, [4]), surface mail communication 
(j8 ps 1.5j5[]), web-browsing (j3 ps 1.2, [6]) and library loans 
(j3 ps 1, [7]). In some other cases a monotonous relation has 
been reported between the user activity and the interevent 
time exponent, for example in short-message communication 
(J3 G (1.5,2.1), [8]) or in watching movies (J3 e (1.5,2.7), 
[91]). In a recent paper there can be found a distribution of ex- 
ponents of various channels of communications ifioll . These 
observations make it necessary to find a model in which the 
interevent time exponent is tunable. 

Similar bursty behaviour has been observed in other natural 
phenomena, for example in neuron-firing sequences ll ill or in 
the timings of earthquakes f 12]. The interevent time distribu- 
tion does not give us any information about the dependency 
among the consecutive actions. Correlation between events 
is usually characterised by the autocorrelation function of the 
timings of detected events. Bursty behaviour is often accom- 
panied by power-law decaying autocorrelation function 111311 
which is usually thought to indicate long-range dependency, 
see e.g. I1 1411 . However, time series with independent power- 
law distributed interevent times show long-range correlations 

linn. 



This letter is organized as follows. We start with introduc- 
ing a task-queuing model which has an advantage compared 
to the Barabasi-model ifnll . namely that the observable is not 
the waiting time of an action (between adding to the list and 
executing it) but the interevent time between similar activi- 
ties. We determine the asymptotic decay of the interevent 
time distribution in a simple limit of the model. We give a 
simple proof of the scaling law between the exponents of the 
interevent time distribution and the autocorrelation function 
based on Tauberian theorems. Finally, we demonstrate that 
the scaling law can be violated if the interevent times are long- 
range dependent. 

The model: We assume that people have a task list of size N 
and they choose their forthcoming activity from this list. The 
list is ordered in the sense that the probability w,- of choos- 
ing the f activity from the list is monotonically decreasing 
as a function of position ;. The chosen activity is going to 
be executed and it jumps to the front of the list (fig[T])- This 
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FIG. 1: Dynamics of the list. In every timestep a random position 
is chosen from the list and the activity sitting at the chosen position 
jumps to the first position and is going to be executed. The other 
activities are shifted to fill in the gap. 

mechanism is responsible for producing the bursty behaviour 
because once a person starts to do an activity, that is going to 
have high priority for a while. We analyse the model in which 
the survival function Q„ = 1 — Y.'k=i w k is power-law decay- 
ing, i.e. Q n = cn~ a+l + 0{nT a ). ff{n m ) means that that term 
is asymptotically at most of the order of n m . In this case w, 
is also power-law decaying with exponent a. Exponentially 
decaying priorities would result in trapped activities, meaning 
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that they get executed many times in sequence. If an activity 
happens to leave the trap it may never be executed again (in 
the infinite system limit). Hence such choice-distribution can- 
not lead to power low decaying interevent time distribution. 

The model is capable of covering many types of activities 
but now we only concentrate on one activity marked with A 
in figQ](e.g. watching movies). For the sake of simplicity we 
assume that the list contains only one item of the activity of 
interest. 

The interevent time is equal to the recurrence time of the 
first element of the list. Numerical results show that the in- 
terevent time distribution is power-law decaying with an ex- 
ponential cutoff. The cutoff is the consequence of reaching the 
end of the list from which a geometrically distributed waiting 
time follows: P w t{t) ~ (1 — Wfif . For the sake of simplicity 
we determine the exponent of power-law decaying region in 
the case of an infinitely long list where the exponential cutoff 
does not appear. 

Solution for N — > °°, ( (7 > I): Let q(n,t) (n>l) denote the 
probability of finding the observed element at position n after t 
timesteps without any recurrences up to time t . The restriction 
not to recur is important because this makes large jumps to the 
front of the list forbidden for the observed element. The initial 
condition is set as q(n, 1) = (1 — wi)8„ 2 and time evolution is 
given by 

9(1,0 = (1) 
q{n,t + \) = (1 -Qn-\)q(n,t) + Q„_iq(n- l,t) (2) 

Our aim is to determine the asymptotic behaviour of 
L~ = i<?(»,0- The main trick we use in analysing the recur- 
sion above is to consider the n = const levels first instead of 
the t = const levels which intuitively one would do. At every 
step the time coordinate gets increased while the position of 
the element might remain unchanged or get increased by one 
as well (fig|2|. The path of the element in the (n,t) plane can 




FIG. 2: Some examples for the path in the (n,t) plane. A typical 
section in the path corresponding to the 4t random variable is em- 
phasized by a rectangle. 

be divided into sections that start with a step on the bias and 
are followed by some steps upwards (this can be zero). These 
sections can be characterised by their height-difference which 
can take values from 1,2, These height differences are in- 
dependent and (optimistic) geometrically distributed with pa- 



rameter depending on the position. Let §j be independent ge- 
ometrically distributed random variables with parameter £>A, 
i.e. P(<^; = t) = (1 - QiY^Qi- With n fixed q{n,t) is the 
probability that we find the element at the n th position at time 
t. This corresponds to paths with n — 1 steps to the right 
(started from position 1 at time 0) and q(n, t ) is the distribution 
mass function of the sum of the random heights. 



q(n,t) 




(3) 



Though the random variables ^ are not identically distributed 
- by somewhat tedious but otherwise standard Fourier ana- 
lytic methods (analysis of characteristic functions) - one can 
show that the central limit theorem holds for this situation. In 
this approximation q(n,t) is Gaussian in variable t with mean 
Y!j=\ and variance YIi=\ ^ 2 (<■>/)■ Using integral approx- 

imation to evaluate the sums and keeping only the highest or- 
der terms in n yields: 
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This formula shows that the probability of finding an element 
at position n at time t is centered on the t = curve which 
is in agreement with the numerical results (fig0). We are in- 
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FIG. 3: Contour plots of q(n,t). The image on the left is a numerical 
result, the right plot shows the CLT approximation. The dashed curve 



terested in the sum of q(n,t) in variable n because that gives 
the survival function of the interevent time distribution. We 
approximate this sum by integral and we apply the following 
substitution (with new variable r): 

n= (crjf) 1 / C7 -r(c(7f) 1/(2<7) 



(5) 



For 7> (here / may refer to: a - 1 /2, (7,2(7- 1): 

X. 2y-l 2y-l 

n' = (cat)" — y(cot)~z°~r+o(t~z°~), 

where o(n m ) means that that term is asymptotically of strictly 
smaller order than n m . From equations (I!]-©) it follows that 



(6) 



s: 



q(n,t)dnfa c ((7f) a a ) 



(7) 
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Differentiating this equation gives the first main result of our 
letter, P ie (t) ~ t~P with J3 = 2 - i 

Another characteristic property of the time series is auto- 
correlation function. LetX(t) denote the indicator variable of 
the observed activity: X(t) = 1 if the observed activity is at the 
first position of the list at time t (i.e. an event happened) and 
X(t) = otherwise. By definition the autocorrelation func- 
tion: 



E[X(0)X(f)]-(E[X(0]>? 



(8) 



The stationary solution of the Markov-chain defined by the 
model is uniform, hence E [X(t)] = and 



(0 = 



p(x(t) = ipr(o) = i) 
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The probability of finding the activity at the first position 
can be calculated numerically by successive application of 
the Markov-chain transition matrix. Numerical computations 
show that the autocorrelation function is power-law decaying 
with an exponential cutoff (figH)). Given <j, the autocorrela- 
tion functions for various list sizes can be rescaled to collapse 
into a single curve (figH] inset). This property can be written 




FIG. 4: Autocorrelation function of the model with various list sizes. 
Inset: the curves corresponding to different Ns collapse into a single 
curve if N S A(t) is plotted as a function of t/N?. 



in a mathematical way: 



/ t 



- ry HNr) (10) 



The exponents used for rescaling the autocorrelation functions 
are listed in tabQ] indicating that in the N — > °° limit the au- 
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TABLE I: Numerical results for data collapse in fig[4] scaling pa- 
rameters y and S for various values of a. 

tocorrelation function is power-law decaying with exponent 



a = ^ for <7 > 2. With the exact j3 = 2 — ^ result this can be 
combined into a scaling law: a + j3 = 2. 

Proof of the scaling law: The essential properties of the 
model for the scaling relation are that the interevent times are 
independent and power-law decaying. Let T denote the set of 
recurrence times and let z be an interevent time. The autocor- 
relation function can be written in the following form: 
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which simplifies to s/{t) = P(f G T) if j8 < 2. The Laplace 
transform of the autocorrelation function can be expressed by 
the Laplace transform of the interevent time distribution: 
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= 1-1 



(12) 



Tauberian and Abelian theorems connect the asymptotics of 
a function with the asymptotics of its Laplace transform lfl9ll . 
Applying Abelian theorem to the right side of the last equation 
results in g(X) ~ A Then applying the Tauberian theorem 
yields st{f) — ?^~ 2 or 



a + p =2 



(13) 



To be precise we had to use an extended version of the Abelian 
theorem which can be derived from the original theorem using 
integration by parts. With similar train of thought the scaling 
law can be extended to the 2 < j3 < 3 region where j3 — a = 2 
holds. These results are in agreement with |3 Till . 

The independence of interevent times was important in the 
proof of the scaling law. As a counterexample we constructed 
a long-range dependent set of interevent times by Metropolis- 
algorithm. The base of the algorithm is constructing a Markov 
chain on the integers that has power-law decaying stationary 
distribution P(x) ~ x~P . The algorithm uses a proposal den- 
sity (transition rate) Q{x' ,x„) which generates a proposal sam- 
ple from the current value of the interevent time. This sample 
is accepted for the next value with probability 



a(x 



Q(x n ,x')P(x') 
Q(x',x„)P(x„) 



(14) 



otherwise the previous value is repeated. To generate a power- 
law distributed sample with long-range dependency the mix- 
ing of the Markov chain should be slow, i.e. the gap in 
the spectrum of the Markov chain should vanish. In this 
order we allow only small differences between consecutive 
interevent times: Q(x l ,x n ) = ^I(x n —D,x n +D) or equiva- 
lently x 1 as a random variable is discrete uniformly distributed: 
x' ~ DU[x„ — D,x„ +D]. The simulation results (figHJl show 
that the autocorrelation function decays slower than it would 
be assumed from the scaling law (0. This means that power- 
law decaying autocorrelation function signs long-range de- 
pendencies only if the exponents violate the scaling relation 
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FIG. 5: Simulation results of the Metropolis algorithm. The exponent 
of the interevent time distribution is set to /3 = 1.5 and parameter D 
of the proposal density is varied. The dashed line shows the decay 
corresponding to the scaling law of independent processes. 



A trivial extension of the model could be putting more items 
of the observed activity into the list. Then this parameter 
would tune the frequency of doing the activity. In this case 
the interevent times become dependent and simulation results 
show that the interevent time distribution decays faster than 
power-law but slower than exponential. It is easy to prove that 
in spite of this the autocorrelation function (0 is independent 
of the number of the observed activity and remains power-law 
decaying. 

Different functions for the choice-probability w, will result 
in different interevent time distributions, however, the method 
we introduced using the central limit theorem is general and 
can be applied in many cases. 

The model may be useful not only for human dynamics but 
also for the theory of card shuffling. We can define the time 
reversed version of the model in which we choose a position 
from the same distribution as in the original model and we 
move the first element of the list to the random position. This 
model has similar properties to the original model, e.g. this 
model has the same interevent time distribution and autocor- 
relation function. If we think of the list as a deck of cards, 
then the time reversed model is a generalisation of the top-in- 
at-random shuffle method. The latter is a primitive model of 
shuffling cards: the top card is removed and inserted into the 
deck at a uniformly distributed random position [20]. 

In a proper model one should take into consideration the 
dependency among the interevent times besides the interevent 
time distribution. In human dynamics the latter is slowly de- 
caying and as a consequence of this the autocorrelation func- 
tion of the interevent times is not well-defined (i.e. the second 
moment does not exist). However, the autocorrelation func- 
tion of the time series exists and it might be a good measure 
for long range dependency among the interevent times. Time 
series of messages sent by an individual in an online commu- 
nity are reported to be not more correlated than an indepen- 
dent process with the same interevent time distribution II 1 611 . 
Similarly, the exponents in the neuron-firing sequence s ap - 
proximately satisfy the scaling law (a = 1.0, /3 = 1.1 IU3IP . 



For these systems the model we studied might be applicable. 
However, this is not always the case. For example, the au- 
tocorrelation function of the e-mail sequence decays slower 
(a = 0.75 [13]) than it should do estimated from the scaling 
law. This indicates long range dependency in the time series 
in addition to the power-law decaying interevent time distri- 
bution. In this case a dependent model should be applied, for 
example a model based on Metropolis algorithm that is sim- 
ilar to the one we studied as a counterexample to the scaling 
law. 

In real networks interactions between individuals have to 
be taken into account to reproduce some social phenomena, 
e.g. temporal motifs [21] observed in a mobile call dataset. 
Interactions could be incorporated in this model by allowing 
the actual activity of an individual to modify the priority list 
of some of his/her neighbour. If this effect is rare and can 
be considered as a perturbation, our results for the dynamics 
of the list could be a starting point for further calculations 
covering for example information flow in a network. 
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